clear all
set more off
global dpath "[Userpath]\data and code"


use "$dpath/data/countyEmission.dta", clear

tostring countycode,g( countycode_str ) 
g citycode_str=substr( countycode_str ,1,4) 
g a_str=substr( countycode_str ,-2,. ) 
destring a_str,g(a) 

gen exempt = 1 if citycode_str=="1521" | 									    ///
					  citycode_str=="2224" | 								    ///
					  citycode_str=="4190" | 								    ///
					  citycode_str=="4521" | 								    ///
					  citycode_str=="4600" | 								    ///
					  citycode_str=="5323" |								    ///
					  citycode_str=="5325" | 							        ///
					  citycode_str=="5328" | 								    ///
					  citycode_str=="5329" | 								    ///
					  citycode_str=="5331" | 								    ///
					  citycode_str=="6541" | 								    ///
					  citycode_str=="6540" |								    ///
					  citycode_str=="6543"
replace a_str="00" if a<20 & exempt != 1
g countycode2_str= citycode_str+ a_str
destring countycode2_str,g( countycode2 ) 
drop countycode_str a_str citycode_str a countycode2_str exempt
order countycode2,after( countycode )


sort countycode2 cic2 year
foreach var in wastewater cod nh3n gas so2 dust {
	bysort countycode2 cic2 year: egen `var's = sum(`var')
	drop `var'
	rename `var's `var'
}
duplicates drop countycode2 cic2 year, force

rename countycode countycode_ini
rename countycode2 countycode
drop dq
order countycode countycode_ini


***************************** match with station info **************************
joinby countycode using "$dpath/data/countyStation.dta", unmatched(both)
tab _merge

keep if _merge == 3
drop _merge

gen citycode = int(countycode/100)
gen provincecode = int(citycode/100)

joinby countycode year using "$dpath/data/countyControls.dta", unmatched(both)
tab _merge

keep if _merge == 3
drop _merge


joinby countycode using "$dpath/data/auto_up.dta", unmatched(both)
tab _merge

drop if _merge == 2
drop _merge
replace autoyear_up = 33000 if autoyear_up == .

gen auto_up = 1 if year>=autoyear_up
replace auto_up = 0 if year<autoyear_up

joinby stationcode using "$dpath/data/auto.dta",unmatched(both)
tab _merge
drop if _merge == 2
drop _merge
g cityboundary1=cityboundary-sea-nationalboundary


****************************** Variable Preparation ****************************
gen llgdpC = log(lgdpC)
drop if llgdpC == .

gen llgdp = log(lgdp)
drop if llgdp == .

foreach var in wastewater cod nh3n gas so2 dust {
	bysort countycode year: egen `var's = sum(`var')
	drop `var'
	rename `var's `var'
	gen l`var' = log(`var')
}
duplicates drop countycode year, force

xtset countycode year

winsor2 lwastewater, replace cuts(0.5 99.5) trim
winsor2 lcod, replace cuts(0.5 99.5) trim
winsor2 lgas, replace cuts(0.5 99.5) trim
winsor2 lso2, replace cuts(0.5 99.5) trim
winsor2 ldust, replace cuts(0.5 99.5) trim

xtset countycode year

gen manualupstream = 1 if manualyear != 33000
replace manualupstream = 0 if manualupstream == .

gen treated = 1 if autoupstream == 1 & year >= autoyear
replace treated = 0 if treated == .

gen treatedM = 1 if manualupstream == 1 & year >= manualyear
replace treatedM = 0 if treatedM == .

cd "$dpath\results"

drop if year>=stationyear
replace treated = 1 if autoupstream == 1 & year >= autoyear-1
g lpop=ln(pop)

sort citycode year
foreach i of varlist wastewater cod pop lgdp{
	by citycode year: egen `i'_city=total(`i')
	g l`i'_city=ln(`i'_city)
}
g llgdpC_city=llgdp_city-lpop_city

reghdfe treated lwastewater lcod llgdp lpop, absorb(countycode year) cluster(citycode)
estimates store reg1
reghdfe treated lwastewater lcod llgdp lpop, absorb(countycode stationcode year) cluster(citycode)
estimates store reg2
reghdfe treated lwastewater lcod llgdp lpop, absorb(countycode year) cluster(countycode)
estimates store reg3
reghdfe treated lwastewater lcod llgdp lpop lwastewater_city lcod_city llgdp_city lpop_city, absorb(countycode year) cluster(citycode)
estimates store reg4
reghdfe treated lwastewater lcod llgdp lpop, absorb(countycode citycode#year) cluster(citycode)
estimates store reg5

esttab reg1 reg2 reg3 reg4 reg5 using TableAIV.rtf,replace b(%9.3f) se scalar(N r2_a) title(Table AIV)  star(* 0.1 **  0.05  *** 0.01)
	   
est clear

